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This article provides an introduction to the ideas behind the multilevel blocking (MLB) 
approach to the fermion sign problem in path-integral Monte Carlo simulations, and also 

^r\ ' gives a detailed discussion of MLB results for quantum dots. MLB can turn the expo- 

nential severity of the sign problem into an algebraic one, thereby enabling numerically 
exact studies of otherwise inaccessible systems. Low-temperature simulation results for 

pH , up to eight strongly correlated electrons in a parabolic 2D quantum dot are presented. 

pH , 1 Introduction: The fermion sign problem 

o. 

O , Quantum Monte Carlo (QMC) techniques are among the most powerful methods 

for the computer simulation of strongly correlated many-fermion systems, capable 
of delivering numerically exact results. This article deals with a finite-temperature 
QMC method, namely path-integral Monte Carlo (PIMC), which is based on a 
discretizcd path-integral formulation of the imaginary-time many-fermion propa- 
gator. Despite its promises, applications of QMC to many-fermion systems have 
been severely handicapped by the infamous 'sign problem'. El Exchange leads to 
non-positive-definite fcrmionic density matrix elements, and the sign cancellations 
arising from sampling fermion paths manifest themselves as a small signal-to-noisc 
ratio, 77 ~ exp(—N/3Eo), that vanishes exponentially with both particle number N 
and inverse temperature [3 = 1/ksT (Eq is a system-dependent energy scale). Apart 
from variational or approximate treatments, such as the fixed-node approximation, 
a solution to the sign .problem in QMC simulations had remained elusive. 
*^J ■ In a recent paper, 13 we developed a general scheme for tackling the sign problem 

in PIMC simulations. The method has been applied to interacting electrons in 
a quantum dot a and to the real-time dynamics of simple few-degrees-of-freedom 
systems. Q This multilevel blocking (MLB) approach is a systematic implementation 
of a simple blocking strategy. The theorem behind the blocking strategy asserts 
that by sampling groups of paths ('blocks') at the same time, the sign problem can 
always be reduced compared to sampling single paths as would be done normally 
(see Sec. EJ for details). By suitably bunching paths together into sufficiently 
small blocks, the sign cancellations among paths within the same block can be 
accounted for without the sign problem, simply because there is no sign problem for 
a sufficiently small system. The MLB approach is then able to turn the exponential 
severity of the sign problem into an algebraic one. In practice, that implies that 
significantly larger systems can be now studied by PIMC. 

The purpose of this article is twofold. First, we want to present the basic ideas 



underlying MLB, omitting technical implementation issues discussed in our original 
work. □ In particular, the reason why blocking paths together is helpful will be 
discussed in some detail. The second aim is to demonstrate the practical usefulness 
of MLB for interesting quantum many-fermion applications. In Sec. [|, MLB results 
are presented for strongly correlated electrons in a 2D quantum dot. 

2 Multilevel blocking approach 

2.1 Blocking strategy 

We consider a many-fermion system whose state is described by a set of quantum 
numbers r denoting, e.g., the positions and spins of all particles. These quantum 
numbers may correspond to electrons living on a lattice or in continuous space. For 
notational simplicity, we focus on calculating the equilibrium expectation value of 
a diagonal operator or correlation function (this can be easily generalized), 

where Y]^t represents either a summation for the case of a discrete system or an 
integration for a continuous system. In PIMC applications, imaginary time is dis- 
cretized into P slices of length e = (3/ P. Inserting complete sets at each slice 
m = 1, . . . , P, and denoting the corresponding configuration on slice m by r m , the 
diagonal elements of the reduced density matrix at r = rp entering Eq. (Fy) read 

p 

r*i,...,rp_i m— 1 

One can then construct accurate analytical approximations for the short-time prop- 
agator. This formulation of the problem excludes effective actions such as those 
arising from an integration over the fermions using the Hubbard-Stratonovich trans- 
formation, El, since they generally lead to long-ranged imaginary-time, interactions. 
The MLB approach suitable for such systems is described elsewhere. H 

In dealing with a many-fermion system, we need to sum over all particle permu- 
tations and the best way to do this is to antisymmetrize the short-time propagators 
explicitly instead of letting the Monte Carlo handle it. This leads to the appearance 
of fermion determinants. Strictly speaking, the antisymmetrization has to be done 
only on one time slice, but the intrinsic sign problem is much better behaved if one 
antisymmetrizes on all time slices. 

Choosing the absolute value of the product of the short-time propagators in 
Eq.(g) as the positive definite MC weight function P[-X"], one then accumulates the 
sign &[X] associated with every path X = (fi, . . . , rp) sampled, 

{) PlTO ' (3) 

Assuming that there are no exclusivity problems in the numerator so that A[X] is 
well-behaved, we can gauge the severity of the sign problem in terms of the variance 



of the denominator, 
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where N s is the number of MC samples taken and the stochastic averages are 
calculated with P as the weight function. For the fermion sign problem, where 
<f> = ±1 and hence (<1> 2 ) = 1, the variance of the signal is controlled by the size of 

Remarkably, one can achieve considerable progress by simply blocking paths 
together, a By this we mean instead of sampling single paths in the MC, we can 
sample sets of paths ( "blocks" ) . Under such a blocking operation, the stochastic 
estimate for (A) takes the form 



(A) 



E B (ExeBP[XMX]A[X]) 
Eb(Ex £B P{XMX}) 



(5) 



where one first sums over the configurations belonging to a block B in a way that 
is not affected by the sign problem, and then stochastically sums over the blocks. 
The summation within a block must therefore be done non-stochastically, or al- 
ternatively the block size must be chosen sufficiently small. Of course, there is 
considerable freedom in how to choose this blocking. 

Let us analyze the variance a' 2 of the denominator of Eq.(|5J). We first define 
new sampling functions in terms of the blocks which are then sampled stochastically, 
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, $'[B]=sgn(E P W <i> [ X ]) 
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Rewriting the average sign in the new representation, i.e., using P'[B] as the weight, 
then inserting the definition of P' and $' in the numerator, 
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and comparing to the average sign in the standard representation using P[X] as the 
weight, we obtain 

\(V)\_ ZxP[X] 

K*)l E b p'[b}' 

By virtue of the Schwarz inequality, 
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we see that for any kind of blocking, the average sign improves (or stays the same) , 
K$') I > I (3?) |- Furthermore, since ($ /2 ) = ($ 2 ) = 1, we conclude from Eq.(0) that 



a' 2 <a 2 



(7) 



and hence the signal-to-noise ratio is always improved upon blocking configurations 
together. Clearly, the worst blocking one could possibly choose would be to group 



the configurations into two blocks, one with positive sign and the other with negative 
sign. In this case, blocking yields no improvement whatsoever, and the '<' becomes 
'=' in Eq.(^). It is apparent from Eq.(|7]) that the blocking strategy provides a 
systematic approach to reduce the sign problem. 

2.2 Multilevel blocking approach 



A direct implementation of the strategy described in Sec. 2.1 does indeed improve 
the sign problem but will not remove its exponential severity. The reason is simply 
that for a sufficiently large system, there will be too many blocks, and once the 
signals coming from these blocks are allowed to interfere, one again runs into the sign 
problem (albeit with a smaller scale Eq). n The resolution to this problem comes from 
the multilevel blocking (MLB) approachB where one applies the blocking strategy in 
a recursive manner to the blocks again. In a sense, we form new blocks containing 
a sufficiently small number of elementary ones, and repeat this process until only 
one block is left. 

Instead of allowing for an uncontrolled interference of the block signals, we 
therefore subdivide the block space by forming a hierarchy of 'levels' I = 0, . . . , L, 
where the Trotter number must be of the form P = 2 L . All elementary blocks 
are distributed onto these levels, and a MC sweep proceeds from the bottom level 
£ = up to the top level £ = L. When sampling on a given level £ < L, no sign 
problem is present if sufficiently small block sizes have been chosen. The nontrivial 
computational task consists of finding a controllable way of transferring interference 
information from level t to t + 1. It is then indeed possible to proceed without 
numerical instabilities from the bottom up to the top level, where the expectation 
values of interest are computed. 

The blocks at the top level still have different signs, but the interference is 
however much weaker than in the original formulation. Empirically, we found that 
the exponential severity is turned into an only algebraic one. The algebraic scaling 
can also be derived by a simple argument, u For lack of space, we refer the inter- 
ested reader to our original papera for practical implementation issues of the MLB 
algorithm, and now turn to an application. 

3 Application: Quantum dots 

Quantum dots are solid-state artificial atoms with tunable properties. Confining a 
small number of electrons N in a 2D electron gas in semiconductor heterostructures, 
novel effects due to the interplay between confinement and the Coulomb interaction 
have been observed experimentally. For small TV, comparison of experiments to 
the generalized Kohn theorem indicates that the confinement potential is parabolic 
and hence quite shallow compared to conventional atoms. Employing the standard 
electron gas parameter r s to quantify the correlation strength, a Fermi-liquid-like 
picture with the sequential filling of single-particle orbitals is applicable only for 
small r s . In the low-density (strong- interaction) limit of large r 8 , however, classical 
considerations suggest a Wigner crystal-like phase with electrons spatially arranged 
in shells, Q termed Wigner molecule due to its finite extent. We are particularly 



interested in the crossover regime between these two limits, where single-particle 
or classical descriptions break down, and basically no other sufficiently accurate 
method is available. Exact diagonalization is limited to very small N since one oth- 
erwise introduces a huge error due to the truncation of the Hilbert space. Hartree- 
Fock (and related) calculations become unreliable for large r s and incorrectly favor 
spin-polarized states. Furthermore, density functional calculations can introduce 
uncontrolled approximations. In fact, by comparing with our numerically exact 
data, symmetry-broken spin density solutions found in a recent density functional 
study, El and later in a Hartree-Fock calculation, E3 arc most likely artefacts of the 
approximations involved. A variational Monte Carlo with a fixed-node approxima- 
tion has been employed by Bolton. O Comparing this to our exact calculations, we 
find that typical fixed-node errors in the total energy for N > 5 are of the order 
of 10%. It is then clear that one must resort to exact methods, especially when 
examining spin-dependent quantities. 

3.1 Model 

We study a clean 2D parabolic quantum dot with zero magnetic field, 
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Here the electron positions (momenta) are given by Xj (pj), their effective mass is 
to*, and the dielectric constant is k. The MLB calculations are carried out at fixed 
N and fixed z-component of the total spin, S = (JVj — N\)/2. We have studied the 
energy, E = (H), the radial charge and spin densities p(r) and s z (r) normalized 
to L dr2nrp(r) — N and J dr2nr s z (r) — S, and the two-particle correlation 
function 

^® = ]vW3T)(£/ (f -^ + ^ ) ) ' (9) 

g s is isotropic, and with y = r/lo prefactors are chosen such that L dyyg s (y) = 1. 
The confinement length scale Iq = yJh/m*u)Q allows the interaction strength to 
be parametrized by A = lo/a = e 2 /ku>oIq, where a is the effective Bohr radius 
of the artificial atom. For any given N and A, the parameter r s = r* /a follows 
by identifying r* with the first maximum in ^2 S g s (r). In the simulations, unless 
noted otherwise, the temperature was set to kuT — 0.1 fuvo- Data were collected 
from several 10 4 samples for each parameter set {N, S, A}, with a typical CPU time 
requirement of a few days (for each set) on a SGI Octane workstation. To check 
our code, we have accurately reproduced the known exact solution for N = 2. 

3.2 Charge and spin densities 

Figure [j] shows the radial charge density p(r) of the spin-polarized state S — N/2 
for N = 5 to 8 electrons. For A = 2, increasing N does not change p(r) qualitatively, 
but the situation is different for strong interactions (A = 8), where one can observe 





Figure 1: Density p(r) of the spin-polarized state (S = N/2) for A = 2 and A = 8. Dashed, solid, 
dashed-dotted, and dotted curves correspond to N = 5, 6, 7, and 8, respectively. Units are such 
that lo = 1. 



shell formation in real space. Such a spatial structure is clear evidence for Wigner 
molecule behavior. The classical shell filling sequence has been computed recently. El 
For N < 6, the electrons arrange on a ring, but the sixth electron then goes into the 
center (shell filling 1-5). Furthermore, electrons 7 and 8 enter the outer ring again. 
These predictions are in accordance with our data. Clear indications of a spatial 
shell structure at N > 6 can be observed even for A = 4, albeit quantum fluctuations 
tend to wash them out. For A >, 4, the charge densities are basically insensitive to 
S. This is characteristic for a classical Wigner crystal, where the Pauli principle 
and spin-dependent properties are of secondary importance. Our numerical results 
for the spin density in this regime simply follow the corresponding charge density 
according to s z (r) ~ (S/N)p(r). A significant S'-dependence of charge and spin 
densities is observed only for weak correlations. 

3.3 Crossover from Fermi liquid to Wigner molecule 

To study the crossover from weak to strong correlations, we employ the 'spin sen- 
sitivity', normalized to unity for r s — 0, 
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The correlation function g s (r) in Eq.(pl) is a very sensitive measure of Fermi statis- 
tics, in particular revealing the spin-dependent correlation hole. As interactions 



tend to wash out the Fermi surface, the quantity (10) is largest for a Fermi gas, 
r s = 0. Since for r s — ► oo, gs(r) becomes completely spin-independent, ^iv(^s) 
decays from unity at r s = down to zero as r s — > oo. The functional dependence 
of this decay provides insight about the crossover phenomenon under study. 

Figure |2| reveals that the function ^Ar(r s ) becomes remarkably universal and 
depends only very weakly on N. Its decay defines a crossover scale r c , where an 
exponential fit for small r s yields r c re 4. For r s > 4, the data can be fitted by 

£(r a ) - exp (- yfcfc) , (11) 



where r' c re 1.2. Remarkably, this is precisely the behavior expected from a semi- 
classical WKB estimate for a Wigner molecule. H The crossover value r c re 4 is also 
consistent with the onset of spatial shell structures in the density, and with the 
spin-dependent ground state energies expected for a Wigner molecule. Therefore 
the crossover from weak to strong correlations is characterized by the surprisingly 
small value r c re 4, instead of r c re 37 found for the bulk 2D electron gas. i3 The 
stabilization of the Wigner molecule can be ascribed to the confinement potential. 
In the thermodynamic limit, uq — > with r s fixed, plasmons govern the low-energy 
physics, and hence the bulk value r c re 37 becomes relevant for very large N. For 
GaAs based quantum dots, we estimate that for 7V,< 10 4 , the va ue r c re 4 is valid. 
Remarkably, very recent experiments on vertical quantum dots 113 have found ev- 
idence for an even smaller crossover scale r c re 1.8. The experimental study was 
carried out in a magnetic field, and the dot contained several impurities. Since both 
effects tend to stabilize a Wigner crystallized phase, our theoretical prediction and 
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Figure 2: Numerical results for £jv(r s ). Statistical errors are of the order of the symbol size. The 
dotted curve, given by exp(— r s /r c ) with r c = 4, is a guide to the eye only. The inset shows tine 
same data on a semi-logarithmic scale as a function of y^rj. The dashed line is given by Eq. (|ll|). 



the experimental observation appear to be consistent with each other, especially 
since we are concerned with a rather smooth crossover phenomenon. 



3-4 Spin- dependent energies 

MLB results for the energy at different parameter sets {N, S, A} are listed in Table 
[D. For given N and A, if the ground state is (partially) spin-polarized with spin S, 
the simulations should yield the same energies for all S' < S. Within the accuracy 
of the calculation, this consistency check is indeed fulfilled. For strong correlations, 
r s > r c , the spin-dependent energy levels differ substantially from a single-particle 
orbital picture. In particular, the ground-state spin S can change and the relative 
energy of higher-spin states becomes much smaller than hug. 

For N = 3 electrons, as r s is increased, a transition occurs from S — 1/2 to 
S = 3/2 at an interaction strength A w 5 corresponding to r s as 8. For N — 4, we 
encounter a Hund's rule case with a small- r s ground state characterized by S = 1. 
^From our data, this standard Hund's rule covers the full range of r s . A similar 
situation arises for N — 5, where the ground state is characterized by S = 1/2 for 
all r s . Turning to N — 6, while one has filled orbitals and hence a zero-spin ground 
state for weak correlations, for A = 8 we find a S = 1 ground state. A similar 
transition from a S — 1/2 state for weak correlations to a partially spin-polarized 
S = 5/2 state is found for N = 7. Finally, for N = 8, as expected from Hund's rule, 
a S = 1 ground state is observed for small r s . However, for A >4, corresponding 



Table 1: MLB data for the energy for various {N, S, A} parameter sets, 
denote statistical errors. 



Bracketed numbers 
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3/2 
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8.37(1) 
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5/2 
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42.86(4) 
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1/2 
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8.16(3) 
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3/2 
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42.82(2) 
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3/2 
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11.05(1) 
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1/2 
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42.77(4) 
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1/2 
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11.05(2) 
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10 


48.79(2) 
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3/2 
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13.43(1) 
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10 


48.78(3) 


3 
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10 
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to r s > 10, the ground state spin changes to S = 2, implying a different 'strong- 
coupling' Hund's rule. 

Let us finally address the issue of 'magic numbers'. For small r s , the filling of 
orbital shells and Hund's rule arguments predict that certain N are exceptionally 
stable. Results for the energy per electron, En/N, in the spin-polarized state 
S = N/2, are shown in Figure 0. Notably, there are no obvious cusps or breaks 
in the ./V-dependence of the energy. Our A = 2 data in Fig. H suggest that an 
explanation of the experimentally observed magic numbers EJ has to involve spin 
and/or magnetic field effects. Remarkably, the absence of pronounced cusps in. 
En/N for strong correlations (A = 8) is in accordance with the classical analysis. 
Therefore magic numbers seem to play only a minor role in the Wigner molecule 
phase. 
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Figure 3: Energy per electron, En/N, for S = N/2 and h(3u>o = 6, in units of hu>o, for A = 2 
(squares) and A = 8 (diamonds). Statistical errors are smaller than the symbol size. Open circles 
are T = fixed-node QMC resultsliil for A = 2. 
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